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1. Introduction 

Inverting the Wilson-Dirac operator is one of the dominant costs in simulating lattice QCD 
with dynamical fermions. A possibility to avoid the explicit inversion is to replace it by a polyno- 
mial approximation of the inverse Wilson-Dirac operator. This allows moreover for a deviation 
from importance sampling with the Boltzmann factor to be compensated by reweighting, which can 
be useful. In addition simulations with odd numbers of flavors become possible. Our focus here 
is the approximate inversion of the non-hermitian Wilson-Dirac operator. According to reference 
[^] the approximation of the non-hermitian operator is supposed to be superior to the hermitian 
version. Each approximation depends on the typical spectrum of the operator to be approximated. 
Therefore, we here report about some spectral properties of the Wilson-Dirac operator with respect 
to Schrodinger Functional(SF)[Q] boundary conditions (BC) and investigate how they affect the 
polynomial approximation. 

2. Computing the spectrum semianalytically 

In the hopping parameter representation Wilson's fermion action ^ reads 

Sf = £ W) [Sxy ~ KHxy] ¥ (y) , (2- 1) 

xy 

and the Wilson-Dirac operator is defined by Mxy = 5^, — KHxy, where the nearest neighbor interac- 
tion is carried by the hopping operator 



Hxy= t (u^)(l-Yfi)8 x+ ^ y + 11^(1 + Yv)8 x - [ i,y). (2.2) 

/i=0 v 



Important spectral features of H depend on the boundary conditions. In general H is not normal 
i.e. [H,W] ^0. 

In order to get a first idea of the spectrum and to develop methods for the general case we start 
by computing the eigenvalues of H in the free case i.e. = 1 in the SF with vanishing background 
field. Since translation invariance in the spatial directions still holds we set 

xir(x) = y(xo)-e m . (2.3) 

using plain waves for the space dependence. This ansatz leads to a reduced 1 -dimensional operator 
of the form 

E = iz»h 0+ i±»ht +i)1O) a 2 = f^sm 2 (p k ) (2.4) 
1 1 k=\ 



acting on y/Q. In (2.4) h is the 1-dimensional hopping operator given for the SF by the nilpotent 
matrix 

(0 1 - 0\ 
-01/ 
- 0/ 
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By itself it has only one (T — 1) fold degenerate zero eigenvalue. The matrix E is related to H by 

3 

spec(//) = 2 spec (E) +2 T\ cos (pk)- (2.6) 

k=\ 

Finding the eigenvalues of E is equivalent to locating the zeros of a smaller determinant 



= det 



A 2 + A (h () +h$) + i4h +<r 



(2.7) 



We could not obtain them in closed form, but approximations are possible. Here we simply com- 
pute the eigenvalues of E numerically for some range of a. 

These eigenvalues are shown in Fig. [T] together with the corresponding ones for (anti)periodic 
boundary conditions that follow trivially from Fourier expansion. For small a the latter approach 
the value 1 leading to zeromodes in 1 — Kji (up to lattice artefacts for antiperiodic BC), where K c 
equals 1/8. In the SF we see that the eigenvalues are 'deflected' away from unity. For very small 
a one can show the behavior 

lAojoca 1 /?- as a ^0. (2.8) 
This is how a gap of order l/T is maintained in the SF (even in the continuum limit). 




Figure 1: Numerical spectrum of E for T = 16 Figure 2: Specttum of the hopping operator// on a 16 4 
and a 2 = tiny ... 3. lattice with = II and SF boundary conditions. 

We conclude this subsection by showing in Fig. ^| a complete spectrum of H computed nu- 
merically. An ellipse with major half axis a = 7.971, minor half axis b = 3.932 and eccentricity 
e = \/« 2 — b 2 = 6.933 is drawn that obviously describes the spectral boundary very well. In the 
centers of the void areas there are (degenerate) eigenvalues for which a vanishes due to the zero 



mode of ho. Roundoff errors in the eigenvalue routine in combination with the behavior (2.8) 
'inflate' three of these dots to small circles. 
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3. Approximation of the inverse Wilson-Dirac operator 

To start simple we first approximate the inverse Wilson-Dirac operator by a geometric series 

M- 1 «P n G (M) = J>»y = i [ _ {k ) h) , for \\kH\\ < 1, (3.1) 
and define the remainder 

R n+1 = II -MP„ G (M) = (KH) n+l . (3.2) 

By construction R n +\ is a small quantity and vanishes in the limit of n — > oo. Moreover, it allows for 
a recursive implementation and its convergence can be easily monitored by computing ||7?, J+ iTj ||, 
where T7 is a Gaussian random vector normalized to 1. Approximating M~ l by a geometric series 
requires a circular bound on the spectrum of radius r = f^An^//)! < 1. From the spectral radius 
r follows the rate of convergence 

V G {K) = -\n{K\X mayi (H)\). (3.3) 

As we have seen in the previous section the shape of the spectrum is elliptical. This fact can 
be exploited to improve our approximation. Expressing the remainder R n+ i in terms of scaled and 
translated Chebyshev polynomials T n [||] we can derive an improved, recursive description, where 
only the eccentricity e as elliptical parameter enters 

R n+l (M) = r "+i((^)/ g ) = a n KHR n (M) + (1 - a n )R n ^ (M) (3.4) 

with R { {M) = kH; R (M) = 1; a n = (l -a„_ie 2 /4) _1 and a x = (l-e 2 /2y\ The second 
equality in ( |3.4[ ) follows from the recurrence relation of the Chebyshev Polynomials, T n+ \(z) = 
2zT n {z) — T„_i(z). By virtue of the defining relation for the remainder we obtain also a recursive 
expression for the Chebyshev approximation of M~ 1 

/f(M) = a n (l + KHP n -\{M)) + (1 -a„)P„_ 2 (M) (3.5) 

with Pi (M) = a\ (1 + kH) and Pq(M) = 1. The rate of convergence in the limit of n — > 00 follows 
using the identity T n (z) = cosh(«arcosh(z)) and by replacing kH by its eigenvalues (cf. footnote 
[[]) [^]. The rate \i c depends on the elliptical parameters a and e which themselves are proportional 
to K 

M C (M=lnf-^p4Y (3-6) 
For periodic BC the extent of the ellipse is known in the free case (a = 8fc, e = \^4Sk) and thus jj, c 



becomes a function of k only, H c (k) = In ^(1 + \[\ -48fC 2 )/(12fc)J , and vanishes like (|3|) for 
K — ► K r . 



'For non-diagonalizable matrices the same behavior is true asymptotically. This can be shown with the help of the 
Schur-decomposition. 
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4. Numerical results 

To test our approximations numerically we monitor the norm of the remainder as a function of 
n and determine the convergence rate from the exponential decay. We perform this test choosing 
n = 400 and varying e to obtain a scan over the eccentricity. In Fig. |3| the results are presented for 
various lattice sizes. 




Figure 3: Rate of convergence pi as function of the Figure 4: Sketch how to get a good guess on the 
eccentricity e. = II, K = 0.115 and each point is eccentricity starting from the largest eigenvalues, 
determined after n — 400 iterations. 

Obviously, there is a dependence on the lattice size in the case of the Schrodinger functional 
and for larger lattices the rate of convergence approaches the value corresponding to periodic BC 
which is independent of L (red dots in Fig. ||). Since the rate of convergence increases roughly by 
a factor 2 from a circular bound (e = 0, geometric series) to the optimal eccentricity, it is important 
to find this optimal value. Therefore we look again at the spectrum of the Wilson-Dirac operator 
and focus our attention especially on the eigenvalue X\ with largest real part and %2 with largest 
imaginary part of H as indicated in Fig. || (multiplied by fc = 0.115). 

One way to obtain a guess on the eccentricity e is to use the norm of A2 as value for the minor 
half axis b. We are then seeking the ellipse which also passes through X\. By the parameter form 
of an ellipse, x = a cos p ; y = b sin p , and using x + iy = X\ we find the major half axis 

a = Re{Ai}/cos(p) with p = arcsin(Im{Ai}/&). (4.1) 

Thus we can determine e = V a 2 — b 2 . Beside yielding a guess on e we can moreover obtain an 
estimate on \i c by eq. (^6) in this way. 



There exist different methods to determine e. In practice we are seeking a good guess on e 
such that the convergence rate is high and its determination is easy. These properties hopefully 
carry over when including a non-trivial gauge field. There we hope to find an eccentricity that 
changes only weakly between different gauge fields at fixed j8 and K. 

For a first experiment with a non-trivial gauge field we start by generating 50 pure-gauge 
configurations on an 8 4 lattice at j8 = 6.0 employing a Cabbibo-Marinari update [§]. Reading these 
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configuration with MAT LAB (version 7.3) and using its implementation of the Arnoldi algorithm we 
try to compute X\ and A2 on each configuration. Unfortunately, the algorithm converged only on a 
subset of the configurations. Hence the mean values presented in Tab. [T] are just a rough estimation 
and within the quoted errors no dependence on the configuration is seen. 





u» 


K 


r 


M G 


a 


b 


e 


M C 


SF 


11 


0.115 


0.8489 


0.1638 


0.9060 


0.4444 


0.7895 


0.1781 


SF 


J8 =6.0 


0.135 


0.838(4) 


0.1770(7) 


0.843(5) 


0.47(1) 


0.701(7) 


0.268(3) 


P 


11 


0.115 


0.9200 


0.0834 


0.9200 


0.4600 


0.7967 


0.1506 


P 


J8 =6.0 


0.135 


0.865(5) 


0.1456(8) 


0.869(6) 


0.48(1) 


0.725(2) 


0.226(9) 



Table 1: Expected values for /I and e derived from measured maximal eigenvalues. SF Schrodinger 
functional, P periodic boundary conditions. 



The results indicate that the spectrum of H for non-trivial gauge fields is expected to be 
"rounder" and enclosed in an elliptical disc of smaller area than the one of the trivial gauge field. 
We check our expectation by computing the polynomial remainder using the above determined e in 
case of the Chebyshev approximation. While for the trivial gauge field we find perfect agreement 
of both methods, the differences in the rate of convergence are larger for non-trivial gauge fields 
when using Chebyshev polynomials 





M G 


e 


M C 


SF 


0.1776(3) 


0.701 


0.2127(2) 


P 


0.1482(4) 


0.725 


0.1914(2) 



Table 2: Computing the convergence from the remain- 
der test as a measure on the "integrated spectrum". 



To get a better understanding of this situation we computed for one SF configuration 800 
eigenvalues with largest/smallest real part and 400 eigenvalues with largest/smallest imaginary 
part again using MAT LAB. Figure |5] shows these data points in blue and the solid red line is the 
ellipse (e = 0.701, a = 0.843) derived from the eigenvalue computation. The predicted value for 
jlc disagrees because the shape of the spectral boundary is not elliptical. The circular bound is 
still estimated correctly as can be seen by the dotted black circle with radius r = 0.838. Hence 
pi G is in agreement. To illustrate our approximation using Chebyshev polynomials we note that by 
specifying e a family of confocal ellipses is determined. From this family the ellipse of smallest 
extent which encloses all eigenvalues specifies a and b, which enter into ([^). Hence the dashed 
red ellipse in Fig. || represents better the one corresponding to the Chebyshev approximation. Here 
a = 0.872 and we compute \i c = 0.208. 

Increasing a to 0.857 and b to 0.547 we yield a different ellipse (e = 0.660) shown with dash- 
dotted green line. This one encloses the computed spectrum even better and leads to the prediction 
\l c {e = 0.660) 0.221. Taking this smaller value of e as input for the Chebyshev approximation 
we find for the rate of convergence \i c = 0.2207(2) confirming the predicted value. 
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e CD 
- e=0.660 
- e=0.701 
e=0.701 



-0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 
Re X 

Figure 5: Computing 2400 eigenvalues of kH on one gauge config- 
uration at j3 = 6.0 and K = 0.135 with SF boundary conditions. 

5. Conclusion and outlook 

These preliminary studies show that a good understanding of the structure of the spectrum 
seems to be important to implement an algorithm approximating the inverse Wilson-Dirac oper- 
ator with good performance. Moreover, useful information on how to tune such an algorithm is 
obtained. 

Probably, the deviation of the spectrum in the SF from an elliptic disk is an artefact of small 
lattice sizes. A check with e.g. a 12 4 lattice would be desirable but seems to be numerically chal- 
lenging. Moreover we like to study the effect of 0(a) improvement (Sheikholeslami-Wohlert term) 
and the effects of preconditioning. 
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